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Spatial Fourier transforms of quasipatterns observed in Faraday wave experiments 
fNl ■ suggest that the patterns are well represented by the sum of 8, 10 or 12 Fourier 

modes with wavevectors equally spaced around a circle. This representation has 
been used many times as the starting point for standard perturbative methods of 
computing the weakly nonlinear dependence of the pattern amplitude on param- 
eters. We show that nonlinear interactions of n such Fourier modes generate new 
modes with wavevectors that approach the original circle no faster than a constant 
times n -2 , and that there are combinations of modes that do achieve this limit. 
As in KAM theory small divisors cause difficulties in the perturbation theory and 
the convergence of the standard method is questionable in spite of the bound on 
the small divisors. We compute steady quasipattern solutions of the cubic Swift- 
Hohenberg equation up to 33 rd order to illustrate the issues in some detail, and argue 
that the standard method does not converge sufficiently rapidly to be regarded as 
a reliable way of calculating properties of quasipatterns. 
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1 Introduction 



Experimental observations of regular patterns have been widely reported in 
many physical systems, for example, Rayleigh-Benard convection, reaction- 
diffusion problems and the Faraday wave experiment [1,2]. In the last example, 
a tray containing a layer of fluid is subjected to vertical vibrations, and the 
flat horizontal surface of the fluid becomes unstable once the amplitude of the 
vibration exceeds a critical value. With multi-frequency forcing, this experi- 
ment is capable of producing a wide variety of patterns with an astonishing 
degree of symmetry [3]. 

The simplest patterns: stripes, squares and hexagons, have reflection, rotation 
and translation symmetries, and a comprehensive theory has been developed 
to analyse the creation of these patterns from the initial flat state [4]. In or- 
der to apply the theory to the experiments, two idealisations are necessary: 
first, the experimental boundaries are ignored, and so in effect the experi- 
ment is supposed to be taking place in a container of infinite size: there are 
two unbounded spatial directions; and second, the observed pattern is sup- 
posed to have perfect spatial periodicity. Restricting to a spatially periodic 
sub domain enables rigorous theory to be applied [5], and the existence of 
stripe, square, and hexagon (and other) solutions of model partial differen- 
tial equations (PDEs) can be proven using equivariant bifurcation theory [4]. 
Given that in some highly controlled experiments the idealisation of spatial 
periodicity appears to hold over dozens of wavelengths of the pattern, these 
assumptions are perfectly reasonable when the objective is to understand the 
nature of these periodic patterns. 

However, experiments are quite capable of producing patterns that cannot be 
analysed in this way. Notable examples of this include quasipatterns, which 
are quasiperiodic in any spatial direction, that is, the amplitude of the pat- 
tern (taken along any direction in the plane) can be regarded as a sum of 
waves with incommensurate spatial frequencies. Experimental photographs of 
quasipatterns are reproduced in figure l(a-c) from [3]; the lack of spatial pe- 
riodicity can be seen in the images, and the long range rotational order is 
evident from the spatial Fourier transforms, in a manner similar to quasicrys- 
tals [6]. The Faraday wave experiment has been a particularly fruitful source 
of quasipatterns since they were discovered by Christiansen et al. [7] (8-fold 
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Fig. 1. (a-c) Experimental observation of (a) 8-fold, (b) 10-fold, and (c) 12-fold 
quasipatterns in the Faraday wave experiment with (a) 3 and (b,c) 2 frequency 
forcing. Top line: experimental photographs; bottom line: spatial Fourier transform. 
The circles in the Fourier spectra indicate wavenumbers that are excited by one of 
the harmonics in the forcing (both harmonics in case b). From Arbell and Fineberg 
(2002) [3], with permission, (d-i) Synthetic examples of periodic patterns (d,e) with 
Q = 4, 6 wavevectors with equal amplitudes, and quasipatterns (f,g,h,i) with Q = 8, 
10, 12, 14 wavevectors. 
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quasipattern with single frequency forcing and a relatively low viscosity fluid) 
and Edwards and Fauve [8,9] (12-fold quasipattern with two frequency forcing 
and a high viscosity fluid). See references [3,10,11,12] for further large aspect 
ratio experiments, and [13] for a recent review of experimental and theoretical 
issues. Quasipatterns have also been reported in studies of nonlinear optical 
systems [14,15], and in numerical studies of several model PDEs [16,17,18], 
though of course computations carried out in large periodic domains can only 
approximate a quasiperiodic pattern. 

The forcing in Faraday wave experiments can be a simple sinusoid, but quasi- 
patterns are more readily generated when two or three commensurate tempo- 
ral frequencies are included in the forcing. Each frequency excites, or nearly 
excites, waves with a particular wavenumber, and nonlinear resonant inter- 
actions between these waves, as well as waves that damped, can be used to 
encourage or discourage particular waves to appear in the pattern [9,19,20,21]. 
By tuning such parameters as the driving frequency, the fluid viscosity and 
the layer depth, experimentalists have been able to produce very clean ex- 
amples of quasipatterns (as in figure la-c), as well as to demonstrate a clear 
understanding of the physical mechanisms behind their production. 

Since quasipatterns, by their very nature, do not fit into periodic domains, 
equivariant bifurcation theory cannot be applied, and other methods are re- 
quired in order to predict, for example, the dependence of the amplitude of 
the quasipattern on parameters, or the stability of the quasipattern. One ap- 
proach, which has been followed many times, is to suppose that the dynamics 
of the quasipattern is dominated by the evolution of the amplitudes of Q waves 
(Q even), with wavevectors distributed equally around a circle. Equations gov- 
erning the evolution of these amplitudes can readily be written down, and take 
the form 



where A, is the complex amplitude of mode j, \i describes the forcing of the 
pattern, and the f3j t k coefficients depend on the angle between the wavevectors 
of modes j and k. Resonant terms (depending on the value of Q) are also 
included but not written explicitly above. Steady solutions of (1) with all 
amplitudes equal represent patterns of the type shown in figure l(d-i), for 
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Q — 4, . . . , 14 wavevectors. 



Equation (1) can be written down directly from symmetry or general physi- 
cal considerations [22,23,24,25], but it has also been derived from PDEs that 
model the hydrodynamic and other problems [26,27,28,29]. The method that 
is used is called modified perturbation theory [30,31], and the amplitude equa- 
tion (1) is only the leading order approximation to the equations that govern 
the evolution of the amplitudes of the modes. 

The problem with this approach when applied to quasipatterns is that it over- 
looks the near-singular resonances that small divisors can cause, and ignores 
the nearly neutral modes that are driven by high-order nonlinear interactions. 

The difficulty of small divisors arises in a variety of situations as well as this 
one, for example, the persistence of quasiperiodic oscillations in Hamiltonian 
systems, the analysis of which cumulated in the KAM theorem (cf. [32]). To 
take an illustrative example, consider the one- dimensional ordinary differential 
equation with quasiperiodic forcing: 

mi,m2=-oo 



-7 



where C mi , m2 are constants satisfying Co,o = and C mijm2 < K± (\mi\ + \m 2 
for every pair of integers m\ and m 2 , with 7 > 2 so that the sum converges. 
The frequencies u)\ and u 2 are incommensurate. This series can be integrated 
formally term by term to give: 

00 /-< 

a U\ — ^mi,m 2 c i(m 1 uj 1 +m 2 uJ2)t ^ 



It is clear that this sum for a(t) may not converge even if the sum for the forcing 
function does, since m\UJ\ + m 2 uj2 comes arbitrarily close to zero, and so the 
amplitudes of the Fourier coefficients can be arbitrarily large. However, if uj\ 
and uo 2 satisfy a Diophantine condition, that is, if there are constants K 2 > 
and 5 > such that cu 1 and u 2 satisfy 

\m1U1 + m 2 co 2 \ > K 2 (|rai| + \m 2 \y 5 (4) 



for every m\ and m 2 , then the sum for a(t) can readily be shown to converge 
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provided 7 > 2 + 5. This increases the constraints on the smoothness of the 
forcing function. Normally, integrating a Fourier series poses no difficulties, but 
this example demonstrates that when quasiperiodic functions are involved, an 
extra degree of caution is necessary. 

Currently, there is no KAM-like theory for quasipatterns, which are quasiperi- 
odic in two dimensions (x, y), rather than than the usual one dimension (time). 
There is, however, a theory for one-dimensional steady quasipatterns [33], 
which makes use of space as a time-like evolution variable. In the absence 
of a rigorous theory for two-dimensional quasipatterns, we examine the issue 
of convergence or otherwise of the standard method of computing amplitude 
equations of the form of (1), when applied to quasipatterns. 

In section 2, we revisit the standard method of modified perturbation theory as 
applied to the computation of the amplitude of a quasipattern as a function of 
a parameter close to the onset of the pattern, and point out where the problem 
of small divisors arises. In section 3, we work out just how small the small 
divisors are (with numerical results that make use of a rapid method presented 
in Appendix A), and in section 4 return to general issue of convergence. We 
discuss the specific example of the Swift-Hohenberg equation in section 5, 
and conclude in section 6 with the argument that the standard method does 
not appear to converge sufficiently rapidly to be regarded as a reliable way of 
calculating properties of quasipatterns. 



2 Perturbation theory 



In order to point out exactly where the difficulties lie, we begin by going 
through the standard modified perturbation theory [30,31] for a general pattern- 
forming PDE: 

^ = ^(U) = C,(U)+^(U), (5) 



where U (x, y, t) represents the order parameter (or any measure of the pat- 
tern), is an operator containing spatial derivatives that depends on a pa- 
rameter ji and that can be split into linear and nonlinear (A/" M ) parts. 
The order parameter may be multi-dimensional. A specific example is the 
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Swift-Hohenberg equation [34]: 

^ = f,U-(l + VrU-U\ (6) 



where U (x, y, t) G R, but many pattern-forming problems can be cast into this 
form, or variations [35]. 

The spatially uniform trivial state U(x, y, t) — is always a possible solution, 
but it loses stability as the parameter \x increases through a critical value /xq, 
which we take to be zero. We focus on the case where the mode that becomes 
unstable is a Fourier mode with nonzero wavenumber, which we scale to 1. The 
problem is posed on the whole plane, so (x,y) G R 2 , but if we were interested 
only in spatially periodic solutions, then the whole plane could be restricted 
to a periodic domain, and standard equivariant bifurcation theory [4] could be 
used. However, this rules out the spatially quasiperiodic patterns of interest 
here. 

Instead of equivariant bifurcation theory, we use the older technique of mod- 
ified perturbation theory [30,31], and suppose that the parameter \i is close 
to its critical value /x = 0, that the amplitude of the solution is small and 
that the pattern evolves slowly. We introduce a small parameter e, scale time 
by e _1 , and write 

U = eU 1 + e 2 U 2 + e 3 U 3 + n = en 1 . (7) 



Then C^U) is of order e, and dU/dT and Af^U) are of order e 2 . In many 
examples, there are additional symmetries in the problem, and it may be 
necessary to scale time by e~ 2 and set fi = e 2 /i2 if it turns out that /ii = 0. 
We focus on the general case here, treating the specific example of the Swift- 
Hohenberg equation in section 5. 

The leading order equation, at order e, is 

£o(tA) = 0. (8) 



In the Swift-Hohenberg example, the operator C would be —(1 + V 2 ) 2 . Since 
fj, = is a bifurcation point, the linear operator C is singular with a circle 
of marginally stable Fourier modes in its kernel: C (e lk ' x ) = whenever k = 
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Fig. 2. (a) Schematic growth (decay) rate of a mode e lk ' x , as a function of k = \k\ 
at ji = 0. Modes with k = 1 are marginally stable, (b) Q = 12 wavevectors on the 
circle k = 1. 

|fc| = 1, (see figure 2a) and so (8) has nontrivial solutions of the form 

^(x,y,0=i;^(0e^-, (9) 



where we can select any Q wavevectors kj, with j = 1, . . . , Q, from the cir- 
cle k = 1 (figure 2b). 

In principle, any set of unit length wavevectors is permitted, though if U is 
required to be real, the negative of each vector must also be included. The 
usual choice for periodic patterns is Q = 2, 4 or 6 equally spaced wavevectors, 
for stripes, squares or hexagons, though Q — 8, 10 and 12 have also been 
used in previous studies and have been observed in experiments. Synthetic 
examples of patterns and quasipatterns with Aj = constant for Q — 4, . . . , 14 
are shown in figure l(d-i). 

At each higher order in e, the equation to solve takes the form 

dU 

C (U n ) = —Af n (Ui, . . . , fn) + (10) 



where the term M n is given by the order e™ part of Taylor expansion of J~n(U) 
in powers of e, so it contains nonlinear terms and the parameter /x x from C^iJJ). 
In principle, the equations can be solved order by order, with each U n deter- 
mined by Ui, . . . , U n -i. At each order, the nonlinear terms A/2, A/3, etc., in- 
volve quadratic, cubic, etc., combinations of the original Fourier modes, which 
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implies that these terms will involve powers of (at most) n in the original am- 
plitudes, and that the Fourier spectrum of Af n contains wavevectors made up 
(at most) of all combinations of up to n of the original wavevectors in the 
set K. If we let 

Q 

k m = ^rrijkj, where m G Z Q , (11) 

3=1 

then the wavevectors in the nonlinear terms at order n are all those with k m 
satisfying \m\ =J2j\ m j\ < n - 

In order to solve the equation (10) at order n, the modes present in Af n are 
divided into two classes. First, if a mode has wavevector on the unit circle, 
then, using the orthogonality in M 2 of Fourier modes with different wavevectors 
(or, more properly, using solutions of the adjoint equation and integrating 
over IR 2 ), the coefficient of this mode on the RHS of (10) must be zero (this 
condition is known as a solvability condition). The reason for this is that 
Co(e lk ' x ) is zero if \k\ = 1, and so such modes are not present in the LHS 
of (10). So, for example, the evolution of the amplitudes Aj(t) is determined 
at second order and takes the general form: 

A j = f j (A 1 ,...,A Q ;» 1 ), (12) 

where the dot stands for evolution on the slow time scale. These evolution 
equations will contain only terms linear and quadratic in the A's. In examples 
with additional symmetry (or with Q ^ 6 and Q ^ 12 equally spaced modes), 
this equation is vacuous, with fx\ = and Ui = 0, and the evolution of the 
amplitudes Aj(t) is determined at third order - see section 5. 

Once the modes on the unit circle have been removed by satisfying a solvability 
condition, all remaining modes in in Af n (making up the second class of modes) 
have wavevectors off the unit circle. For these modes, C (e lk ' x ) is nonzero, so 
the singular linear operator operator C can be inverted to give U n : 

U n = Co 1 (-K(U U . . . , C/ n _ i; ^) + ) . (13) 

Inverting the operator C generates arbitrary linear combinations of modes in 
its kernel that are used to satisfy solvability conditions at higher order. Like 
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the nonlinear term Af n , each U n will include modes with wavevectors made up 
of (at most) of all combinations k m of up to n of the original wavevectors, 
with \m\ < n. 

As well as the method outlined above, there are two other approaches to these 
computations, both of which avoid adding modes on the unit circle at each 
order when inverting Cq. First, the original paper [30] used an expansion for 
the parameter \i: \i = e/ii + e 2 /i 2 + . . ., and chose values of /ii, /i 2 , etc., to 
satisfy the solvability conditions at each order. As a result, for a given value 
of /x, a polynomial for e must be solved before the amplitude can be computed. 
In the second alternative, the original \i = e/ii is used, leaving modes in the 
RHS of (10) on the unit circle, which cannot be removed. The modes involved 
will be exactly the Q modes that were taken in the original ansatz for U\, 
and the contributions that appear in (10) can be redesignated as order e n_1 
corrections to the leading order solvability condition (12). We have checked 
for some specific examples that the three approaches give the same results, 
but prefer the approach described in detail for the problem at hand. 

In many cases, the leading order solvability condition (12) is sufficient, but 
in other problems, this equation is degenerate, and the calculation is carried 
to some higher order N in powers of e. An implicit assumption is that the 
power series expansion (7) for U converges for some nonzero e as this process 
is repeated (and N — > oo). At each order n in the perturbation calculation, 
the operator Cq must be inverted for each mode e tkm ' x , where \m\ < n and 
k m could be close to the unit circle. In typical pattern forming problems, the 
growth rate of a mode e lk ' x has a quadratic maximum at k = 1 (see figure 2a), 
so 



for \k\ close to 1, with equality in the case of the Swift-Hohenberg example. 
For periodic patterns (with Q = 2, 4 or 6 modes), integer combinations of 
the initial wavevectors form a lattice, so the wavevectors k m cannot come 
arbitrarily close to the unit circle (apart from the modes on the unit circle, 
which are dealt with by applying solvability conditions). In this case, conver- 
gence will not be a problem for small enough e. However, for quasipatterns 
(with Q — 8, 10, 12 or more modes), there is no lattice and combinations of 
modes can come arbitrarily close to the unit circle, leading to small divisors 




(14) 
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in the denominator when £ is inverted. The issue of convergence has never 
been properly examined for two-dimensional quasipatterns, and most authors 
assume that the leading order solvability condition (12) yields useful and re- 
liable information about the amplitude and stability of the quasipattern. We 
will see below that it is far from obvious that this is the case. 

This so-called small divisor problem is well known in other situations that 
feature quasiperiodicity, and, in particular, is known to arise with these quasi- 
patterns. What is not known is how rapidly Cq 1 grows as the order of trunca- 
tion N increases. The question here is how close k m can get to the unit circle 
as \m\ = N becomes large, and does the power series expansion for U have a 
nonzero radius of convergence in spite of the bad behaviour of Cq 1 . We turn 
to these two questions in the next sections. 

3 Combinations of modes 

In this section, we take integer combinations of up to N of the Q original 
vectors on the unit circle, and compute how close these combinations can get 
to the unit circle as N becomes large. We are able to prove that the closest 
\k m \ can get to 1, with |m| = N, is bounded above and below by a constant 
times N~ 2 in the cases Q — 8, 10 and 12. We also have numerical evidence 
that the closest distance is bounded above and below by a constant times iV -4 
in the case Q — 14 (at least for N < 1000), going to zero much more rapidly 
than in the other three cases. 

We begin with figure 3(a-f), illustrating the locations of combinations of up 
to iV = 7, 11 and 15 wavevectors in the case Q = 12. Note how the density 
of points increases with N, and how the minimum distance between points 
and the unit circle goes down with N. Figure 3(g-i) compares with the cases 
Q = 8 and 10, which show similar behaviour, and with Q = 14, which has a 
much higher density of possible wavevectors for the same value of N. 

Figure 4 shows detailed numerical results for the smallest nonzero distance 
||fe m | — 1| from the unit circle as a function of the total number of modes 
\m\ = N for Q = 8, 10, 12 and 14 original modes. The calculations are made 
possible by a rapid method of searching for the closest approach, presented 
in appendix A: the method is order iV 2 for Q — 8, 10, 12 and order N 4 for 
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(a) 



(b) 







Fig. 3. Positions of combinations of up to N wavevectors original vectors on the unit 
circle, with (a) Q = 12, AT = 7, (b) Q = 12, AT = 11, (c) Q = 12, N = 15; (d,e,f) on 
second row: details of first row. The circle indicates the unit circle, k = 1, the large 
dots are the original Q wavevectors, and the small dots are integer combinations of 
these, (g), (h) and (i) show N = 15 and Q = 8, 10 and 14. Note how the density 
of points increases with iV and with Q, and the proximity of points to the unit 
circle decreases with N. Note also how the density of points is markedly higher 
with Q = 14, for the same value of N. 

Q = 14. The solid lines in figure 4 confirm numerically that the scaling for the 
distance to the unit circle is order iV -2 for Q — 8, 10 and 12, and order iV -4 
for Q = 14. The remainder of this section is devoted to proving the correctness 
of the iV -2 scalings, and in particular, to showing how for certain values of N, 
wavevectors close to the unit circle can be found explicitly in the cases Q — 8, 
10 and 12, using continued fraction expansions. 

We label the vectors k\, k 2 , . . . , kg anticlockwise around the circle start- 
ing with ki = (1,0), with kj + Q/ 2 = —kj. We are interested in the scal- 
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10° 10 1 10 2 10 3 10 4 10 ; 

N 




N 



Fig. 4. Smallest nonzero distances from the unit circle ||fc m | — 1| as a function of 
the total number of modes |m| = TV, for (a) Q = 8, (b) Q = 10, (c) Q = 12 
and (d) Q = 14. Stars in (a-c) mark distances calculated from equations (24-26), 
and straight lines indicate the scaling iV~ 2 . The two staircase-shaped lines in (a) 
indicate minimum distances for even and odd values of N. The straight lines in (d) 
indicate iV -4 . 

ing behaviour of how close k m = J2 ( j=i m jkj can get to the unit circle as 
| m.j = J2j \ m j\ — N becomes large, so we seek the vector of integers m with 
\m\ = N that yields the vector k m that is closest to the unit circle for this 
value of N. Once we have found a particular vector that is close to the unit 
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circle, we are also interested in finding the smallest N that can achieve this dis- 
tance. Including equal and opposite vectors kj and kj + Q/ 2 will only increase N 
without coming any closer to the unit circle, so we take only mi, ... , mq/ 2 but 
allow these to be negative. 

With this restriction, the squared length of a vector k m is, for each value of Q: 



Q 


= 2 : 


\k I 2 


= m 2 


Q 


= 4 : 


\k I 2 


2 i 2 

= m x + m 2 


Q 


= 6 : 


\k I 2 


= ml + ml + m 2 + mim 2 + m 2 m 3 — 7773777 1 


r\ 


o . 

— o . 


\h I 2 


= 777 2 + ml + 777 2 + 777 4 

+ \/2(^i^2 + rn2 m 3 + rn 3 m 4 — m 4 mi) 


Q- 


= 10 : 


\k I 2 


= ml + ml + m 2 + m\ + m 2 

+ ijj{mim 2 + m 2 m 3 + m 3 m 4 + m 4 m 5 — m 5 mi) 

+ (c<j — 1) (777.1777.3 + m 2 m 4 + m 3 m 5 — m 4 mi — m 5 m 2 ) 


Q- 


= 12 : 


\k I 2 


— m\-\- m\ + 7773 + 777,4 + m 5 + m 6 

+ 777x7773 + ?™2 m 4 + m 3 m 5 + ^4^6 ~ m 5 mi — 777 6 777 2 

+ ^(777 1 7772 + m 2 m 3 + m 3 m 4 + m 4 m 5 + m 5 m & — m§mi) 


Q- 


= 14 : 


\k I 2 

| "'ml 


= ml + ml+ 777 2 + 7774 + 777 2 + 777g + 777 2 

+ cji(777i777 2 + m 2 m 3 + 777 3 777 4 + m 4 m 5 + m 5 m e + m e m 7 — 

+ CJ2(^1^3 + 777 2 777 4 + 77737775 + 777 4 777 6 + 777 5 7777 — 7776777 1 — 

+ ^3(777x7774 + rn2 m 5 + m 3 m 6 + m 4 m 7 — m 5 mi — m 6 m 2 — 



777 7 777i) 
77777772) 

m 7 m 3 ) 



where a; is the golden ratio: = (1 + \/5)/2 = 2cos(7r/5), with tu 2 = + 1, 
and o;-,- = 2 cos(j7r/7), with c^2 = c^ 2 — 2 and cj 3 = 1 — u;i +w 2 - The irrational u;i 
is the root of a cubic equation: u;^ — uj\ — 2uj x + 1=0. 

We observe that for Q = 2, 4 and 6, \k m \ 2 is an integer, so, as expected, 
points on square and hexagonal lattices cannot come arbitrarily close to the 
unit circle without actually hitting it. 

For Q — 8, 10 and 12, \k m \ 2 is of the form: 

\k m \ 2 = 1 +p-rq, (15) 

where r = a/2, u or \/3 is an irrational root of a quadratic equation with 
integer coefficients, and p and q are integers, lip — rq is close to zero (that is, 
if r is well approximated by |), then \k m \ 2 can come close to 1. The particular 
rational approximations involved are: 
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I = 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


r = V2 


EL — I 

qi i 


3 
2 


7 

5 


17 
12 


41 
29 


99 
70 


239 
169 


577 
408 


1393 
985 


3363 
2378 


8119 
5741 


r = ijj 


pi — i 

Ql 1 


2 
1 


3 
2 


5 
3 


8 
5 


13 

8 


21 
13 


34 
21 


55 
34 


89 
55 


144 

89 


r = s/3 


EL — I 


2 
1 


5 
3 


7 
4 


19 
11 


26 
15 


71 
41 


97 
56 


265 
153 


362 
209 


989 
571 



Table 1 

Continued fraction approximations to r = y/2, u> and \/3, as a function of the order I 
of the truncation. 



g = 8 

Q= 10 
Q = 12 



Q 



m% 



1 



UJ PS - = 



m 4 mi — m 3 m 4 — m 2 m 3 — mim 2 
m\ + . . . + m 5 m2 — 1 



m^mi — ... — mim 2 



- m 6 m 2 - 1 



m 6 mi 



(16) 
(17) 
(18) 



In the expressions above, we choose p and g, which depend on the integers m, 
to be positive and to have no common factors. 

For Q — 14 (and higher), \k m \ 2 involves the sum of an integer plus integers 
times at least two different irrationals that will not, in general, be roots of 
quadratic equations with integer coefficients. Apart from the numerical evi- 
dence in figure 4, we do not pursue the cases Q > 14 further here. 

It is clear that the theory of continued fraction approximations of irrationals 
will be useful here. The continued fraction expressions for the irrationals r = 
y/2, oo and \/3 are (respectively): 



1 + 



1 + 



1 + 



(19) 



2 + 



2 + 



1 + 



1 + 



1 + 



2 + 



1 + 



2 + 



If these fractions are truncated after I terms, the successive fractions pi/qi that 
approximate r are given in table 1. We recall from the theory of continued 
fractions [36] that 



Ki 
if 



< 



— r 



< 



K2 
if 



(20) 



15 



for r = y/2, lj, s/3, and K 1: K 2 are constants. The values of K x and K 2 , 
which are order unity, are related to the largest integers appearing in the 
expansions (19), which are 1 or 2 in these cases. These inequalities mean that 
the truncated continued fraction expansions ^ approximate r well, but not 
too well, as I becomes large. It will also be useful to note that if / > 1 and if 
q is an integer with < q < qi, then [36] 



; 2 i) 



Pl 


< 


P 


r 


r 


Qi 




Q 



This means that if — is the truncation of the continued fraction approximation 
of an irrational r, no other fraction with a smaller denominator comes closer 
to r. 



If we exclude those vectors k m that fall exactly on the unit circle, which would 
have p = q = 0, the relations in (20) and (21) can be used to show that \k m \ 2 
can approach 1 no faster than order N~ 2 for Q — 8, 10 and 12. The reason 
is that the denominators in (16-18) are of the form of a sum of products 
of the integers m 1? . . . , toq/ 2 . Since \m,j\ = N, each quadratic term in the 
denominator can be no larger than N 2 in magnitude; there are no more than Q 
of these terms, so the denominator as a whole satisfies q < QN 2 . Then, using 
first (21) and then (20), we have 



k r - 1 



\p — rq\ > \pi — rqi\ > 



Ki 
Qi 



(22) 



where qi is the smallest of the g/'s above q: < q < qi (unless — qi — 1). 
Now for r = v^2, w and a/3, the g/'s are no further than a factor of 3 apart [36], 
so qi < 3(ft_i < 3q < 3QN 2 , and we have (assuming |fe m | ^ 1): 



> 



K 

N 2 ' 



(23) 



where \m\ — N and K is a constant. 

We now show that the order iV~ 2 rate of approach is indeed achieved for 
Q = 8, 10 and 12. Observe that, for Q = 8: 



km = fci +pik 3 + q t k 6 + qikg = (l,p t - y/2q t ), 



(24) 
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with \m\ = N = pi + 2qi + 1 (even) and \k m \ 2 — 1 = pf + 2qf — 2p l q t y/2 (a 
related vector with N = 2pi + 2qi + 1 (odd) can also be found); for Q = 10: 



k m = (pi + l)fe 2 + (pi ~ l)k 5 + (qi + l)fc 8 + (q t - l)fcg (25) 
= (1, (pj - wgj)V3- w), 

with |m| — N — 2pt + 2q t and |fc m | 2 -l = 3p 2 + 2g 2 + 2piq t - (pf -qf +Ap l q l )ou, 
using the fact that cu 2 — uj + 1; and for Q = 12: 

fcrrr = Pi&4 + (ft ~ 1)^9 + (ft + 1)^11 = - V^ft), (26) 



with \m\ = N = pi + 2q\ and (fc^l 2 — 1 = pf + 3qf — 2piqi\f?>. These vectors 
were found after a prolonged examination of the distances plotted in figure 4. 

Using (20), we have, for these particular vectors, 

= (pi-r qi y<!$, (27) 




where r stands for \/2, to or v^, with an extra factor of 3 — cu in the case r = cu. 
Using relations like qi < pi < ft+i < 3qi [36], and the relations between N, qi 
and pi above, it is then possible to show in each case that iV is less than a 
constant times qi, so 



k r - 1 



K' 



where K' is a constant. These particular choices of k m are plotted on the 
graphs in figure 4(a-c) as stars. 

From the graphs it is clear that these particular vectors are not always the 
closest ones that can be found for given values of iV (particularly for Q = 10), 
but they suffice to prove the scaling results required here. If \k m \ 2 — 1 as given 
above (that is, \k m \ 2 — 1 = pf + 2qf — 2piqi\^2 for Q — 8, and so on) can 
be written as py — rqy for some integer I', then we would expect \k m \ 2 to be 
particularly close to 1, compared with other vectors of up to that order. So, 
for example, an excursion into numerology suggests the following relations, 
which we have proven for Q = 8 by induction: for Q — 8, I' — 21 + 1; and for 
Q = 12, 1' = 21 + 1 if I is odd. On the other hand, for Q = 10, pf - qf + Ap^ 
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does not appear to equal qy for values of / up to 15, which is probably why the 
constructed vector does not achieve the closest possible distance in this case. 



The summary result of this section is that we have shown that, given an 
integer N, the vector k m with \m\ = N that comes closest to the unit circle 
(without being on the unit circle) satisfies 



K 



< 



< 



(29) 



for constants K and K', for Q — 8, 10 and 12 equally spaced original vectors. 
The numerical evidence, for N < 10 6 , suggests values for Q = 8: K = 0.72 and 
K' = 16.95; for Q = 10: K = 0.69 and K' = 9.94; and for Q = 12: K = 0.56 
and K' = 4.34. 

It should be emphasised that several of the steps in the derivation of these 
bounds (for example, equation (20)) rely on the fact that \/2, oo and y/3 are 
quadratic irrationals, that is, they are roots of quadratic equations with integer 
coefficients. This implies that the integers in the continued fraction expansion 
of these numbers (19) form repeating sequences and so are bounded above. 
The case of Q — 14 (and higher) is more difficult to analyse because the 
irrational numbers U\ and 0J2 in the expression for \k m \ 2 in this case are not 
quadratic irrationals, and because there are two irrationals. 



4 The question of convergence 



The results of the previous section, taken with (14), imply that when k m 
is close to the unit circle, C,Q 1 (e lkm ' x ) can be as large as a constant times 
N 4 e ik m X for Q = 8, 10 and 12, with N — \m\. 

At each stage in the calculation, the linear operator C is inverted as in (13): 
U n = Co 1 (-K{U 1 ,...,U n - 1 )), (30) 



where we simplify the discussion by dropping the parameter and the time 
dependence, and we assume that the solvability condition has been satisfied 
so that £ can be inverted . The nonlinear terms Af n will contain modes e tkm ' x 
with \m\ up to and including n, so, at least at first glance, U n will be n 4 times 



18 



larger than the product of various combinations of U\, . . . , U n -\- In particular, 
the combination C/iC/n-i occurs in M n . This suggests that U n is n 4 times larger 
than C/„_i, which is (n — l) 4 times larger than £/ n _ 2 , and so on - so U n will 
be of order (n!) 4 . The perturbation method should yield the pattern U(x,y) 
as the limit: 



but if U n is as large as (n!) 4 , the limit will not converge for nonzero e. For the 
limit to converge, U n should be no larger than a constant raised to the power 
of n. 

In fact, the (n!) 4 estimate is unduly pessimistic, since there will be cancel- 
lations. Moreover, if we focus on nonlinear interactions in J\f n that result in 
modes with wavevectors on the unit circle, the modes that are close to the 
unit circle involved in these nonlinear interactions will originate from U n /2, not 
from U n -\. The reason is that the terms in Af n involve products like £/i?7 n _i, 
U<iU n -2i ■ ■ ■ , Now the term U\U n -i cannot generate modes closer to the 
unit circle than the modes already in U n -\. Modes that are closer than any 
previous combination will come (if at all) from combinations like U^, 2 , and 
will involve the modes in U n /2 that are close to the unit circle. This leads to 
an estimate of the form: U n will be n 4 times larger than U^> 2 , which will be 
(n/2) 8 times larger than £/ 4 / 4 , and so on. The resulting estimate for U n will 
not be as large as (n!) 4 , but it will still be larger than any constant raised to 
the power of n. 

Unfortunately, it is difficult to be more precise than this, but the arguments 
above certainly cast doubt on the convergence of the modified perturbation 
theory method as applied in this way to determine the amplitude of a quasi- 
pattern as a function of parameter. At best the series will be asymptotic, that 
is, not converging and yet still yielding useful information when truncated at 
low order. 

As a cautionary tale, we turn briefly to the asymptotic series representation 
of the Stieltjes integral, as described in [37]: 



N 



U= lim Ye n a 



(31) 




OO oo 



(32) 
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/ is a perfectly well defined integral, depending on a small positive parameter e. 
The fraction may be expanded formally as a power series in e, but this step 
is invalid, as the sum does not converge if t > 1/e, and the upper limit of 
the integral is t — infinity. The formal power series can be integrated term by 
term, resulting in an infinite sum for the value of / as a function of e. This 
sum does not converge for any nonzero e, and yet, if the sum is truncated at 
a fixed order N, then there is a range of e for which the truncated sum gives 
a reasonable approximation of I [37]. Taking a larger N results in a smaller 
range of e, closer to zero. On the other hand, if e is fixed at a small number, 
then there is a a truncation N that gives an estimate of I that is reasonably 
close to the correct value. The truncation can be increased as e is taken to be 
smaller. 

The danger is that the relation between the usable range of e and the level 
of truncation is not known in advance. The most severe truncation (I — 1) is 
the safest, but loses useful information about the dependence of / on e. 

In the next section, we take the specific example of steady quasipattern solu- 
tions of the cubic Swift-Hohenberg equation (6). The equivalent severe trun- 
cation would have the amplitude of the quasipattern be V //I, which has been 
widely used by authors who are relying on the small divisor problem alluded 
to above not rendering this truncation meaningless. 



5 Example: the Swift— Hohenberg equation 

In this section, we go through the details of deriving expressions for the am- 
plitudes in the specific example of steady solutions of the Swift-Hohenberg 
equation (6). This is one of the simplest pattern- forming PDEs, and serves to 
illustrate the problem at hand. 

For this presentation, we concentrate on the case Q = 2, but we have carried 
out the computations for Q — 2, ... ,12 and up to 33 rd order in e. The symme- 
try U — > —U implies that all even terms U2, U4, etc., are absent, that = 0, 
and that time should be scaled by <r 2 . To simplify the presentation, we will 
only seek steady solutions, and so drop the time derivative terms. By taking 
e 2 to be the bifurcation parameter, we can set 112 — 1 (taking /i 2 to be positive 



20 



since the bifurcation is supercritical). The expansion is then 
U = eU! + e 3 U 3 + e 5 f/ 5 + • • • , fi = e 2 , 



(33) 



with C (U) = -{l + V 2 ) 2 U. 

The leading order equation at order e is Cq{U\) — 0, which is solved by (for 
real U and with Q = 2 modes, for clarity of exposition): 

U 1 (x,y) = A 1 e lx + A 1 e- lx , (34) 



where altering the phase of the complex amplitude A 1 translates the pattern. 
At third order in e, the equation we have to solve is: 

C (U 3 ) = -U x + U 3 . (35) 



Now Uf = A\e 3lx +3\Ai\ 2 Aie lx +c.c. (where c.c. stands for complex conjugate), 
and since £0(^3) cannot contain e lx and e~ tx , equation (35) gives: 

= A 1 -3\A 1 \ 2 A 1 . (36) 



This equation would usually have a factor of /x 2 in front of the linear term, but 
we have set this to I. The nontrivial solution of this equation is A 1 = l/y/3, 
where we may take A 1 to be real without loss of generality. The operator £ 
acts on a mode e lk ' x as C (e tk ' x ) = —(1 — k 2 ) 2 e lk ' x , so £0(^3) can be inverted 
to give: 

U 3 = X —= (e 3ix + e" 3 ") + A 3 e lx + A 3 e~ ix , (37) 

192y3 

where we have included an arbitrary combination of the original modes on the 
unit circle (in the kernel of £ ), which will allow the solvability condition to 
be satisfied at the next order. 

At fifth order in e, the equation we have to solve is: 



C (U 5 ) = -U 3 + 3U 2 U 3 (38) 
= —^-^ ( elX + ^ + eblX ) + + ( A 3 + A 3 )e lx + c.c. 
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Table 2 

Values of the coefficients A n in (39). These data are also plotted in figure 5 

The solvability condition can be satisfied by setting A 3 = l/384\/3, with 
an arbitrary imaginary component (set to zero) that corresponds to a small 
translation of the original pattern. With the e lx component removed, £ can 
be inverted to give U$. 

This procedure is repeated up to some order N, resulting in a power series for 
the original pattern U : 
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Fig. 5. (a) Absolute values of the coefficients A n in (39) for Q = 2 (lowest line), 
4, 6, 8, 10 and 12 (top line) modes. These data are also given in table 2. Straight 
lines indicate that the coefficient of the modes on the unit circle will converge as 
the order of truncation iV is increased, and the inverse of the slope of the straight 
line gives the radius of convergence. The dotted lines are 0.3™ (lowest), l n and 4.8™ 
(highest), (b) Ratios A n+ 2/A n against n: this ratio continues to increase with n for 
Q > 8, indicating that the sum in (40) may not converge for nonzero e. 
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where A n is the coefficient of (e/v^3) n in the expression for the amplitude of 
the modes on the unit circle. The factor is chosen so that A 1 — 1 (so 
the Ai and A3 here are rescaled from the A\ and A3 above). The calculation 
can be carried out for other values of Q; for Q — 4, 6, . . . , 12, the scaling for 
amplitudes is \/9, \/l5, . . . , \/33, so A x = 1 in all cases. The values of the first 
few coefficients A n are given in table 2 and plotted in figure 5. 

In these calculations, all modes generated by nonlinear interactions were kept 
for Q = 2, 4 and 6. For the other three cases, we kept only modes with 
wavenumbers up to v^5, to keep the total number of modes within manageable 
limits. Even so, with Q = 12, there were more than 15000 modes generated 
at the highest order - without this truncation, there would have been almost 
2 million. We checked the lower orders against calculations keeping all modes, 
and the differences in the mode amplitudes were about 1%. Restricting the 
number of modes in this way had no effect on how close combinations of 
wavevectors could get to the unit circle. 

For Q — 2, 4 and 6, it appears that there will be no problem with the con- 
vergence of the coefficient of the modes on the unit circle in (39) since the 
coefficients A n grow no faster than a constant to the power of n (indicated 
by straight dotted lines in figure 5a). The bound on the ratios A n+ 2/A n also 
indicates the range of values of e for which convergence is expected. 

For Q — 8, 10 and 12, the situation is less clear. The values of the coefficients 
in table 2 are certainly very large (10 50 ), but they are multiplied by the small 
parameter e raised to high powers. At first glance, the values of the coefficients 
in figure 5(a) appear to be going up as straight lines (which would indicate 
convergence), but the ratios A n+ 2/A n in figure 5(b) show that the slopes of 
these lines plateaus for n in the range 12-18, and then start to increase grad- 
ually for larger n, suggesting that the coefficients A n are growing faster than 
a constant raised to the power of n, and casting doubt on the convergence 
of (39) as the level of truncation is increased. 
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The fact that the ratios A n+2 /A n plateau and then start climbing in figure 5(b) 
can be related to the nonlinear interaction of modes at each order in the 
perturbation theory. Taking Q = 8 to be specific, for 7 < n < 9 and for 
11 < n < 17, the values of \k m \ 2 generated up to these orders, closest to the 
unit circle, are 0.85786 and 1.05887 respectively. Even though these modes 
close to the unit circle are generated at some order n, they do not influence the 
value of the coefficient of modes on the unit circle until nonlinear interactions 
work their way back: this will occur at order 2n + 1. As a result, we expect the 
jump in k 2 from n — 9 to 11 to influence A 23 - indeed, the ratios start to climb 
at A 23 , and more sharply at A 25 . We would then expect these ratios to plateau 
and then start climbing at order 39, plateau and climb again as wavevectors 
generated at higher order come closer to the unit circle. The same reasoning 
would predict the ratios to start climbing at n = 27 for Q = 10 and at n = 15 
and 23 for Q = 12, at least roughly consistent with the data in figure 5(b). It 
was only by going to such high order that these issues became clear. 

Equation (39) can be used to find the amplitude A^ of modes on the unit 
circle as a function of the original parameter jj, = e 2 when the expression is 
truncated to include powers in e up to and including e^: 

AW = y/ji £ A n ^ 2 , (40) 

ra=l,3,... 

where the amplitude has been rescaled (as indicated above) so that A^ = ^fji 
for all values of Q. Graphs of A^ as functions of /i for iV = 1 to 31 and for 
Q = 2 to 12 are presented in figures 6 and 7. For Q = 2, 4 and 6 (figure 6a- 
c and figure 7a), the curves of A^ against /i converge as iV increases in a 
manner consistent with the straight line increases of A n in figure 5(a). In fact, 
an estimate of the radius of convergence can be made from the slopes: the 
series will converge for /i < 32 for Q = 2, fx < 9 for Q = 4, and fi < 0.65 for 
Q = 6. These limits are roughly half the values of fi at which modes generated 
in nonlinear interactions become linearly unstable. 

On the other hand, for Q > 8 (figure 6d-f and figure 7b-d), at each level of 
truncation N, the graph of A^ against /i diverges at a value of /i that is 
decreases as iV becomes larger, consistent with the steady increase in the ra- 
tio A n+2 /A n as n increases. Indeed, modes generated in nonlinear interactions 
become linearly unstable for /i arbitrarily close to zero, for large enough N. 
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Fig. 6. Scaled amplitude as a function of \i from (40), for (a) Q = 2, (b) Q = 4, 
(c) Q = 6, (d) Q = 8, (e) Q = 10, (f) Q = 12, for different levels of truncation 
N = 1, . . . , 31. For Q = 2 and Q = 4, increasing the order of truncation has little 
effect for /x up to 1, while for Q = 6, it appears that increasing the order of truncation 
converges to a solution only for /x < 0.65 or so. However, for Q > 8, increasing the 
order of truncation leads to graphs of as a function of that appear to diverge 
for ii closer and closer to zero as N becomes larger. Details of (c-f) are shown in 
figure 7. 
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Fig. 7. Detail of figure 6, for (a) Q = 6, (b) Q = 8, (c) Q = 10, (d) Q = 12. Note 
how for Q > 8, there is no sign that the graphs of against \i are settling down 
as N increases. 

However, it appears from these graphs that for \x sufficiently small (say, less 
than 0.01 for Q — 8), there is a chance of convergence. Without proper esti- 
mates of the rate of increase of A n with n, it is impossible to know for certain 
whether or not this method converges, and if it does not, whether or not the 
truncated estimates converge to the true solution in the limit of small e. 



6 Discussion 



In summary, we have shown that modes generated by nonlinear interactions 
between Q — 8, 10 and 12 Fourier modes with wavevectors equally spaced 
around the unit circle have wavevectors that can approach the unit circle no 
faster than a constant times n" 2 , where n is the number of modes involved. 
We have also shown by construction that there are combinations of modes 
that do achieve this limit. 
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When carrying out modified perturbation theory in order to compute the 
amplitude of a pattern as a function of a parameter, the usual approach is to 
start with two assumptions: first, that when the parameter is small, the desired 
pattern U can be written as a power series in that small parameter; and 
second, the primary modes of interest have wavevectors equally distributed 
around the unit circle. At each order n in the theory, nonlinear terms generate 
modes involving up to n of the Q modes. The modes that fall exactly on the 
unit circle are dealt with by applying a solvability condition, while equations 
for modes off the unit circle are satisfied by inverting the linear operator £, . 
In the cases Q = 2, 4 and 6, the patterns are spatially periodic and modes 
generated by nonlinear interactions do not approach the unit circle. For Q — 8, 
10 and 12, the wavevector can come within n~ 2 of the unit circle, and small 
divisors (order n -4 ) appear when inverting Co, leading to numerically large 
coefficients in front of the Fourier modes. These coefficients grow sufficiently 
rapidly with n that convergence of the power series for the pattern U is called 
into question. 

We have explicitly carried out modified perturbation theory up to 33 rd order 
for the cubic Swift-Hohenberg equation. Of course, this kind of calculation 
cannot demonstrate convergence or otherwise, but it does illustrate the issues 
that arise. The main conclusion of the calculation is that even if modified 
perturbation theory does generate a convergent series approximation to the 
quasipattern for small enough /x, the series certainly diverges if the parameter \x 
is bigger than about 0.01, depending on exactly which value of Q is used. It 
is possible that the series do converge for smaller /i, though we have argued 
that this is not the case. Even if the series do diverge for all nonzero /i, a low- 
order truncation may still give a useful approximation of the quasipattern, 
assuming that the equations do have a quasipattern solution. It is on this 
basis that other researchers have proceeded. 

There are two related issues at stake. First, existence: do pattern forming 
PDEs like the 2-dimensional Swift-Hohenberg equation have quasipattern so- 
lutions that bifurcate from the trivial solution? Second, given the small divisor 
problem, can asymptotic methods like modified perturbation theory yield use- 
ful approximations to these solutions? We have not addressed the first issue 
in this paper, but plan to turn to it in future. The limits we have derived 
on the rate of approach of wavevectors to the unit circle will play a central 
role in that calculation. As for the second issue, we have shown that modi- 
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fied perturbation theory does not converge sufficiently rapidly (or slowly) to 
provide an answer unequivocally one way or the other, and so this standard 
method should not be regarded as a reliable way of computing properties of 
quasipatterns. 

What is needed is a method that converges more rapidly. Each order in the 
standard theory gains a factor of e as well as large factors from any small 
divisors that arise. There are other methods, developed for proofs of KAM 
theory (see [32]) that converge more rapidly, and these may be required for a 
rigorous treatment of quasipatterns as well. The difference between the KAM 
situation and that of quasipatterns is that in the KAM case, the solutions of 
interest are quasiperiodic in only one dimension (time), while in the second, 
quasipatterns are quasiperiodic in two space directions. 

By making arbitrarily small perturbations to the Q wavevectors, it is possible 
to make qualitative alterations to the nature of the problem in the cases Q — 8, 
10 and 12. For instance, the patterns can be made periodic on square or hexag- 
onal lattices, with a lower limit to how close vectors can get to the unit circle. 

2 2 

For example, in the case Q = 12, choosing the modes (1, 0), ( JPffi , p>~^ ) and 
so on. where ^ is a continued fraction approximation to v^3 (see table 1) yields 
12 modes on the unit circle that become nearly equally spaced as I increases, 
and that generate a square lattice by virtue of (pf — qf,2piqi,pf + qf) be- 
ing Pythagorean triplets - see [38] for more details. Similarly, 12-dimensional 
representations of the group D e ix T 2 can be chosen so that the modes are 
nearly equally spaced and yet they generate a hexagonal lattice [19]. Even in 
a square periodic domain, approximate Q = 12 quasipatterns can be gener- 
ated [24] . The 8-dimensional representations of D4 ix T 2 can be used to ap- 
proximate 8-fold quasipatterns in the same ways, though it is not clear how a 
10-fold quasipattern could be approximated by a periodic pattern. The draw- 
back with approximating quasipatterns by periodic patterns in these ways is 
that the range of validity of the normal forms derived shrinks to zero as the 
approximation improves. 

It is interesting to note that 8, 10 and 12-fold quasipatterns have been observed 
experimentally for several years now, but no 14-fold (or higher) quasipattern 
has been reported (cf. [23]), with the possible exception of [14]. One might 
speculate that the reason for this is that the convergence issues discussed 
above are likely to be more serious in the case Q — 14 since wavevectors 
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approach the unit circle much more rapidly than in the cases 8, 10 and 12 (see 
figure 4). 
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A Appendix: method for finding the closest modes 

In this appendix, we present an order N 2 algorithm for finding which combi- 
nations of up to N vectors end up near the unit circle. The method is suit- 
able for Q — 8, 10 and 12, and can be extended to an order N 4 method 
for Q = 14. We focus on Q = 12 for definiteness, and let k\ = (1,0), 
k 2 = (cos(27r/12),sin(27r/12)), etc. 

For each value of N, we want to find non- negative integers rrij (j — 1, . . . , 12) 
such that k m = J2j m j^j is close to the unit circle, with \m\ = J2j \ m j\ — N, 
and m achieves this minimum distance to the unit circle for this N. In fact, 
we are interested in all combinations with \m\ < N satisfying this minimality 
condition. The requirement of minimality and the symmetries of the problem 
lead to restrictions on the integers rrij that allow the order N 2 algorithm 

By rotating the vectors, we may choose m± > 0, without loss of generality. 

The requirement for minimality amounts to considering only those rrij where 
there is no set of values m'j such that J2j m jkj = J2j m 'jkj and < 
\ m j\- Using mod 12 arithmetic, this leads to: 

• If rrij > then rrij+e = 0, since otherwise let m'j = rrij — rrij + Q and m'j +6 = 
(or the other way round if rrij+Q > rrij), with m[ = mi for I ^ j , j + 6. This 
gives a smaller set of vectors summing to the same point. 

• If rrij > then rrij + 4 = 0, since kj + kj + 4 = kj + 2, and we can let: 
m 'j+2 = m j+2 + min(mj, mj +4 ) m'j = mj — min(mj,mj +4 ) m'j +A = mj +4 — 
min(TOj, mj +4 ) m\ = m; for I ^ j, j + 2, j + 4. 

• Similarly, if mj > then mj +8 = 0. 

Combined with m x > 0, these imply m 5 = m 7 = m 9 = 0. Also, only one of 
m 3 and m u can be nonzero, by the same arguments. Using mirror symmetry, 
choose m 3 > and tou = 0. 

Similar arguments applied to m-i, m^, m$, m§, m±o and m±2 imply that only 
two of these can be nonzero, and those two must be separated by 2: r«2 and 
m 4 , or m 4 and m 6 , or m 6 and m 8 , or m 8 and m w , or m w and m 12 , or m 12 and 
m 2 . Consider these in turn, with mi > and m 3 > 0. 

• m 2 > and m 4 > 0: if these two are nonzero, then all possible values of 
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k m = J2j m jkj are in the upper right quadrant and cannot be close to the 
unit circle. 

• m 4 > and m 6 > 0: all possible values of k m lie in the upper half-plane, 
and cannot be close to the unit circle. 

• mio > and myi > 0, or myi > and m 2 > 0: all possible values of k m lie 
in the right half-plane, and cannot be close to the unit circle. 

The only remaining possibilities are either m 6 > and m 8 > 0, or m 8 > 
and m w > 0. However, for every combination of rrij in one configuration 
(mi, m 3 , 777.6, m s), there is an equivalent combination in the other configuration 
(mi, m 3 , m 8 , m 10 ) that has equal distance to the unit circle: rn'i = m 3 , m 3 = 
mi, m 6 = mio, m' 8 = m 8 , m' 10 = m 6 . So we need only consider cases where 
(mi, m 3 , m 8 , m 10 ) are nonzero. 

Looping over all possible combinations of (mi, m 3 , m 8 , m 10 ) with mi + m 3 + 
m 8 + mio < ^ gives an order iV 4 algorithm, but this can be improved as 
follows. 

The vectors ki = (1, 0) and fci = (0, —1). If k m is to be close to the unit circle, 
m 3 k 3 +m 8 k 8 must lie in or near the upper left quadrant of the wavevector plane 
(positive k y , negative k x ), or at least within the range k x < 2 and k y > —2. 
Furthermore, for given m 3 and m 8 , the values of mi and m w for which k m can 
lie close to the unit circle are quite restricted: mifci + mi fci = (mi, —mio) 
must be near the vector — m 3 fc 3 — m 8 fc 8 . 

So instead of looping over all possible combinations of (m 1 ,m 3 ,m s ,m 1Q ), it is 
only necessary to loop over (m 3 ,m 8 ) and check values of mi close to (within 
2 of) the negative of the x component of m 3 fc 3 + m 8 fc 8 , and values of mi 
close to (within 2 of) the y component of m 3 k 3 + m 8 k s , which results in an 
order iV 2 algorithm. 

The algorithm can be tidied up a little, and similar arguments can be applied 
in the cases Q = 8 and (with a little more difficulty) Q = 10. When Q > 14, 
only order iV 4 or slower algorithms are possible, based on the same ideas. 
These methods were used to generate the data in figure 4. 
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